//----------------------------------*-C++-*----------------------------------//
/*!
 * \file   testMoments_Field.cc
 * \author Jeremy Roberts
 * \date   Jul 11, 2011
 * \brief  Test of moments field.
 * \note   Copyright (C) 2011 Jeremy Roberts. 
 */
//---------------------------------------------------------------------------//
// $Rev::                                               $:Rev of last commit
// $Author::                                            $:Author of last commit
// $Date::                                              $:Date of last commit
//---------------------------------------------------------------------------//

#include <iostream>
#include <iomanip>
#include <cstdio>
#include <cstdlib>
#include <cmath>

#include "test/Testing.hh"

#include "tran/Moments_Field.hh"
#include "tran/CartMesh.hh"
#include "tran/Traits.hh"
#include "tran/Definitions.hh"
#include "../utilities/Constants.hh"
#include "../utilities/DBC.hh"
#include "../utilities/Soft_Equivalence.hh"
#include "../utilities/SP.hh"
// Equations
#include "tran/DD_1D_Equations.hh"
#include "tran/DD_2D_Equations.hh"
#include "tran/DD_3D_Equations.hh"



using namespace std;
using namespace slabtran;
using namespace util;
using namespace constants;


int test1d_1();
int test2d_1();
int test3d_1();
int test1d_index();
int test2d_index();
int test3d_index();

int main(int argc, char *argv[])
{
  Require(argc==2);
  int test = atoi(argv[1]);
  cout << " ------ testMoments_Field ------ " << endl;
  int test_value = 0;
  int number_of_tests = 7;
  if (test == 1 || test == number_of_tests) test_value += test1d_1();
  if (test == 2 || test == number_of_tests) test_value += test2d_1();
  if (test == 3 || test == number_of_tests) test_value += test3d_1();
  if (test == 4 || test == number_of_tests) test_value += test1d_index();
  if (test == 5 || test == number_of_tests) test_value += test2d_index();
  if (test == 6 || test == number_of_tests) test_value += test3d_index();
  if (test > number_of_tests || test < 1) cout << " test " << test
      << " not available. " << endl;
  if (test_value) cout << " failed with value = " << test_value << endl;
  else
    cout << " passed." << endl;
  return test_value;
  return 0;
}

/*!
 *  This tests the generation of a moment's index for all
 *  dimensions.  For each case, the moments are cycled
 *  through by their cardinal index, e.g for 3d
 *  \verbatim
 *    cardinal:   0       1      2     3      ...
 *       (l,m): [(0,0)],[(1,-1),(1,0),(1,1)], ...
 *  \endverbatim
 */
int test1d_index()
{
  int legendre_order = 3;
  // 1d
  for (int l = 0; l <= legendre_order; l++)
    TEST(l == Moments<_1D>::index(l));
  cout << " ... passed 1-d moments index test. " << endl;
  return 0;
}

int test2d_index()
{
  int legendre_order = 3;
  int k = 0;
  for (int l = 0; l <= legendre_order; l++)
    for (int m = -l; m <= l; m++)
      if ((l + m) % 2 == 0)
      {
        TEST(k == Moments<_2D>::index(l, m));
        k++;
      }
  cout << " ... passed 2-d moments index test. " << endl;
  return 0;
}

int test3d_index()
{
  int legendre_order = 3;
  int k = 0;
  for (int l = 0; l <= legendre_order; l++)
    for (int m = -l; m <= l; m++, k++)
      TEST(k == Moments<_3D>::index(l, m));
  cout << " ... passed 3-d moments index test. " << endl;
  return 0;
}

int test1d_1()
{

  // Make a simple mesh.
  double width = 10.0;
  int count = 10;
  Mesh<_1D>::SP_mesh mesh;
  mesh = new Mesh<_1D>(width, count);

  // Initialize the moments field with 0.
  int legendre_order = 2;
  Moments_Field<DD_1D_Equations,_1D>::SP_moments_field phi;
  phi = new Moments_Field<DD_1D_Equations,_1D>(mesh, legendre_order, 0.0);

  // Set to 10*i + m*m = (0,1,4),(10,11,14),...
  // Rememember, moments are cell->unknown->moment
  int k = 0;
  for (int i = 0; i < count; i++)
    for (int l = 0; l < legendre_order + 1; l++, k++)
      (*phi)[k] += 10.0 * i + 1.0 * l * l;

  TEST(k == phi->size());

  // verify
  k = 0;
  for (int i = 0; i < count; i++)
  {
    int cell = mesh->convert(i);
    for (int l = 0; l < legendre_order + 1; l++, k++)
    {
      int moment = Moments<_1D>::index(l);
      double phi_ref = 10.0 * i + 1.0 * l * l;
      // Test operator(cell,unknown,moment) interface.
      TEST((*phi)(cell, 0, moment) == phi_ref);
      // Test operator(cell,unknown,l,m) interface.
      TEST((*phi)(cell, 0, l, 0) == phi_ref);
      // Test moment(i,u,l) interface.
      TEST(phi->moment(i, 0, l) == phi_ref);
    }
  }

  cout << " ... passed 1d moments field test. " << endl;
  return 0;
}

int test2d_1()
{

  // make a simple mesh
  double width = 10.0;
  int count = 10;
  Mesh<_2D>::SP_mesh mesh;
  mesh = new Mesh<_2D>(width, count);

  // initialize the moments field with 0
  int legendre_order = 2;
  int number_moments = 6; // (2+1)(2+2)/2 = 6
  // total size should be (10*10)*6 = 600
  Moments_Field<DD_2D_Equations,_2D> phi(mesh, legendre_order, 0.0);

  /*!
   * The following is the suggested way to update the moments
   * within the equation solver (e.g. \ref DD_2D_Equations).
   * The programmer implementing a 2d solver is \em obligated
   * to update only those moments valid for 2d, which means
   * \em not trying to update \f$\phi^l_0\f$ for \f$l>0\f$. This
   * condition is easily established if the 2d moments interface
   * is used, as it requires the caller uses valid indices; hence,
   * the programmer will be quickly alerted to any errors.  So,
   * use the 2d interface!
   *
   */
  int k = 0;
  for (int j = 0; j < count; j++)
  {
    for (int i = 0; i < count; i++)
    {
      for (int l = 0; l <= legendre_order; l++)
      {
        // for m from -l to l by 2
        for (int m = -l; m <= l; m += 2)
        {
          phi[k++] += (10.0 * i + 100.0 * j - 1000 * m + 10000.0 * l);
        }
      }
    }
  }

  TEST(k == phi.size());
  TEST(legendre_order == phi.legendre_order());
  TEST(number_moments == phi.number_moments());

  // verify
  for (int i = 0; i < count; i++)
  {
    for (int j = 0; j < count; j++)
    {
      int cell = mesh->convert(i, j);
      for (int l = 0; l <= legendre_order; l++)
      {
        for (int m = -l; m < l; m += 2)
        {
          int moment = Moments<_2D>::index(l, m);
          double phi_ref = 10.0 * i + 100.0 * j - 1000 * m + 10000.0 * l;
          // Test operator(cell,unknown,moment) interface.
          TEST(phi(cell, 0, moment) == phi_ref);
          // Test operator(cell,unknown,l,m) interface.
          TEST(phi(cell, 0, l, m) == phi_ref);
          // Test moment(i,j,u,l,m) interface.
          TEST(phi.moment(i, j, 0, l, m) == phi_ref);
        }
      }
    }
  }

  cout << " ... passed 2d moments field test. " << endl;
  return 0;
}

int test3d_1()
{

  // make a simple mesh
  double width = 10.0;
  int count = 10;
  Base_Mesh::SP_mesh mesh; // use base pointer
  mesh = new Mesh<_3D>(width, count);

  // initialize the moments field with 0
  int legendre_order = 2;
  int number_moments = 9; // (2+1)^2 - 2 = 9 - 2 = 7
  // total size should be (10*10)*7 = 700
  Moments_Field<DD_3D_Equations,_3D> phi(mesh, legendre_order, 0.0);

  int n = 0;

  for (int k = 0; k < count; k++)
  {
    for (int j = 0; j < count; j++)
    {
      for (int i = 0; i < count; i++)
      {
        for (int l = 0; l <= legendre_order; l++)
        {
          for (int m = -l; m <= l; m++)
          {
            double tmp = 1.0;
            if (m!=0) tmp = m/fabs(m);
            phi[n++] += tmp * (10.0 * i + 100.0 * j + 1000 * k
                + 10000 * fabs(m) + 100000.0 * l);
          }
        }
      }
    }
  }

  TEST(n == phi.size());
  TEST(legendre_order == phi.legendre_order());
  TEST(number_moments == phi.number_moments());

  // verify
  for (int i = 0; i < count; i++)
  {
    for (int j = 0; j < count; j++)
    {
      for (int k = 0; k < count; k++)
      {
        int cell = mesh->convert(i, j, k);
        for (int l = 0; l <= legendre_order; l++)
        {
          for (int m = -l; m <= l; m++)
          {
            int moment = Moments<_3D>::index(l, m);
            double sign = 1.0;
            if (m != 0) sign = m / fabs(m);
            double phi_ref = sign * (10.0 * i + 100.0 * j + 1000 * k + 10000
                             * fabs(m) + 100000.0 * l);
            // Test operator(cell,unknown,moment) interface.
            TEST(phi(cell, 0, moment) == phi_ref);
            // Test operator(cell,unknown,l,m) interface.
            TEST(phi(cell, 0, l, m) == phi_ref);
            // Test moment(i,j,u,l,m) interface.
            TEST(phi.moment(i, j, k, 0, l, m) == phi_ref);
          }
        }
      }
    }
  }

  cout << " ... passed 3d moments field test. " << endl;
  return 0;

}

//---------------------------------------------------------------------------//
//              end of testMoments_Field.cc
//---------------------------------------------------------------------------//
